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Abstract. We propose a novel approach based on a Langevin equation for fluctuating motion of the center of mass of granular 
media fluidized by energy injection from a bottom plate. In this framework, the analytical solution of the Langevin equation 
is used to derive analytic expressions for several macroscopic quantities and the power spectrum for the center of mass. In 
order to test our theory, we performed event-driven molecular dynamics simulations for one- and two-dimensional systems. 
Energy is injected from a vibrating bottom plate in the one-dimensional case and from a thermal wall at the bottom in the 
two-dimensional case. We found that the theoretical predictions are in good agreement with the results of those simulations 
under the assumption that the fluctuation-dissipation relation holds in the case of nearly elastic collisions between particles. 
However, as the inelasticity of the interparticle collisions increases, the power spectrum for the center of mass obtained by 
the simulations gradually deviates from the prediction of theoretical curve. Connection between this deviation and violation 
of the fluctuation-dissipation relation is discussed. 
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INTRODUCTION 

Granular materials fluidized by external driving have been widely studied in recent years in the field of nonequilibrium 
statistical physics. In order to fluidize granular matter, we have to supply sufficiently large energy continuously, 
because the kinetic energy of grains is dissipated due to inelastic collisions among grains. There are various ways 
to supply energy to the system: a way commonly used in experiments is use of a vibrating bottom plate. In computer 
simulations, a stochastic thermal wall is also often used as an ideal model of energy supply from boundaries. When 
the balance between energy input and dissipation is achieved, the system is in a nonequilibrium steady state. 

It is known that the granular fluids under external driving show various fascinating phenomena, such as convection 
inside the fluids, waves and patterns that appear on the surface, size segregation like the Brazil nut phenomenon, a 
transition from a condensed state to a fluidized state (See Ref. [1] and references therein). It is also important to 
study how macroscopic quantities in the system, e.g., height of the center of mass (COM) and granular temperature, 
depend on the system parameters, such as the number of particles, the restitution coefficient, the amplitude and 
frequency of the external vibration. Experimental studies have revealed some scaling relationships for the macroscopic 
quantities [2, 3, 4, 5]. 

There have been also a lot of theoretical studies to understand these phenomena and properties. At the most 
microscopic scale, molecular dynamics (MD) simulations that follow the motion of all particles using computer 
have been commonly used [2, 3, 6, 7, 8, 9, 10, 11]. One of the most important tools to study dilute granular gas is 
kinetic theory, which is formulated in terms of the Boltzmann (-Enskog) equation under the assumption of "molecular- 
chaos". It has been used to study velocity distributions and the scaling relationships in granular fluids under external 
vibration [4, 12, 13, 14, 15, 16]. At more macroscopic scale, hydrodynamic equations which are derived from the 
Boltzmann equation have been applied to driven granular systems to study local density, flow, and temperature 
profiles [17, 18, 19,20, 21]. 

These theoretical approaches have successfully described a variety of phenomena mentioned above, within the range 
of validity of each approach. It was found, however, some discrepancies in the results of these approaches, when they 
were applied to the scaling relationships for the macroscopic quantities. For example, as summarized in the review 
of previous works by McNamara and Luding [9], experimental results by Warr, Huntley, and Jacques [4], suggest a 
scaling law for the height of the COM h c . m . measured from its height at rest h c . m .o taking the form 

h c .m. -h c . m . ~ {A Q co) s H v , (1) 
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with 8 = 1.3 ± 0.04 and v = 0.27 ±0.11, where A Q a> represents the maximum velocity of the vibrating bottom, and 
H is the number of layers of particles at rest. Results of MD simulations of a two-dimensional system by Luding, 
Herrmann, and Blumen [6], however, suggest 8 ~ 1.5 and v ~ 1 for simulations without rotation of particles; results 
by Luding [7] give 8 — 1.60 ±0.10 and v = 0.76 ±0.11 for simulations including rotation of particles. On the other 
hand, results of kinetic theoretical approaches [4, 13, 14] give 8 = 2, v = 1. Although there are some studies that 
take into account an effect of the wall [9], an effect of high density [15], and an effect of nonuniformity of the 
hydrodynamic fields [20], this discrepancy has not yet been fully explained. 

Recently we have proposed a novel theoretical approach [22, 23] which describes the system at more macroscopic 
level than the previous approaches mentioned above. Our theory is on the basis of a Langevin equation of motion of the 
COM, which is derived by focusing on the force exerted to the granular media by the bottom plate. Some macroscopic 
quantities can be easily obtained from the solution of this equation. In this proceedings, we will first describe how 
to formulate the Langevin equation for the COM, and then test the theory by comparing some important predictions 
derived from the Langevin equation with the result of extensive event-driven MD simulations. A detailed discussion 
on validity of the fluctuation-dissipation relation which connects the intensity of the random force with the friction 
coefficient will be one of the main points of this article. 



THEORETICAL FORMULATION 

In the following we explain our theoretical formulation in the case of D-dimensional system of grains on a vibrating 
bottom plate. The case of grains on a heat source at bottom can be discussed only by a slight modification as shown 
later. 

We consider granular medium that consists of N inelastic particles of mass m and diameter d bouncing on a vibrating 
bottom plate, subjected to gravity with acceleration g. The z-direction is chosen to be opposite to the direction of 
gravity. The motion of particles is confined in a box of cross section S, where the quantity S is an area for D = 3, 
a length for D = 2, and unity for D = 1 . The bottom plate oscillates sinusoidally with amplitude Aq and angular 
frequency CO. Hence, the height of the bottom plate zo{t) at time t is 

z Q (t) =A sin((Of). (2) 

In this paper, we focus only on a fluidized state without any structure in the horizontal directions, such as convection 
and surface wave. Such a fluidized state can be achieved assuming sufficiently small length scales of the cross section 
of the box containing granular medium, so that convections and surface waves are suppressed. Moreover, we ignore any 
boundary effects associated with the side walls for simplicity. Indeed, in event-driven MD simulations to be described 
later, periodic boundary conditions in the horizontal directions have been used. 



Langevin equation 

We will begin by considering the equation of motion of the COM of the particles: 

d 2 Z 

M-^- = -Mg + F b . (3) 

Here, M = Nm is the total mass of particles and Z is the height of the COM. Ignoring the boundary effects concerning 
the side walls, the force acting on the COM is the sum of gravity and the force exerted by the bottom plate F b . 

There are two important timescales in the system: the oscillation period of the bottom plate t(= 2n/a>) and the 
macroscopic relaxation time x re i to the stationary state. If % is comparable to x re i, the energy supplied by one stroke 
of the bottom plate is almost dissipated during the period, and the particles will be in a condensed state. Such a 
condensed state is beyond the scope of the present study. In this study, we restrict ourselves to the high-frequency case 
(T/T re i) 2 <C 1, in which the system is in a fluidized state. 

In order to evaluate F b , the force exerted on the particles by the bottom plate, we consider its reaction force, that 
is, the force exerted on the bottom plate by the particles. Here we draw an analogy between our system and a system 
of a Brownian particle exhibiting Brownian motion (see, for example, Ref. [24]). In our system, granular particles 
randomly collide with the bottom wall and exert a force on it. In Brownian motion, solvent molecules randomly 
collide with a Brownian particle and exert a force on it. These two forces are expected to have similar properties. The 
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force on one side of the Brownian particle consists of a force due to the static pressure, the frictional force, and the 
random force. Thus, on the basis of the analogy, it is plausible to assume that F\, consists of the following four kinds 
of forces: average force due to the static pressure, frictional force, elastic force, and random force. 

The elastic force has its origin in elastic oscillation excited by the bottom plate. An elastic oscillation mode that 
shows the slowest relaxation is the mode with the largest-wavelength in the system. In the stationary state, this 
mode plays a dominant role in deciding macroscopic property of the system. In the largest-wavelength mode, the 
granular fluid shows macroscopic oscillation alternating between the expansion state with the highest surface of the 
fluid and the contraction state with the lowest surface, which accompanies oscillation of the height of the COM with 
the same frequency. It must be noted that the oscillation of the bottom plate does not trigger resonances with the 
largest-wavelength mode and the other macroscopic modes, because the condition we assumed above, (xjx re i) 2 <C 1, 
suggests that the time scale of the bottom plate oscillation is much shorter than time scales of macroscopic motions. 
Another important elastic oscillation is a sound wave that is directly excited by the bottom plate. Hence, we take into 
account two forces as the elastic force: a force due to this largest-wavelength mode and a force due to a sound wave 
excited by the bottom plate oscillation. 

Finally, we assume the following form of Fb(t): 

F b (t)=Mg-MnV z (t) -Mil 2 (Z(0 -Z) +f s +R(t). (4) 

The first term is the average force acting on the bottom plate. The second term is the frictional force, which is assumed 
to be linear in the z-component of the velocity V z of the COM, with a constant coefficient jx . The third term is the elastic 
force resulting from the largest-wavelength mode showing macroscopic expansion and contraction of the granular 
fluid. We assume that it is linear in the height of the center of mass Z measured from its long time average Z. A 
constant coefficient £2 specifies the angular frequency of macroscopic oscillation of the granular fluid. The fourth term 
is the other elastic force resulting from excitation of a sound wave by the vibrating bottom plate. We will formulate 
this term below. The last term is the random force. Its property will also be specified later. 

The period of the macroscopic slowest oscillation T osc is given by z osc = 2n/Cl. Similarly, the macroscopic relaxation 
time x re i is given by T re / = \j\l. Since the both time scales characterize macroscopic changes that extend to the whole 
system, it is plausible to assume that they are on the same order as the characteristic time for a sound wave to travel 
along the vertical direction from the bottom to the surface of the granular gas. Let us define the thermal velocity c, 
by c = yjDksT jm, where T is the global temperature related to the mean square velocity fluctuation and 1cb is the 
Boltzmann constant. Then, the characteristic time can be estimated as c/g because the velocity of sound is on the order 
of c and the height of the surface of the granular gas is, as the first order approximation, on the order of c 2 /2g, which 
is the maximum height of a particle launched from the bottom with the thermal velocity c. Thus, we assume that £2 
and jj. are on the order of g/c: 

£2 = £2^, Ai = £ £ , (5) 

c c 

where £2 and (X are numerical factors on the order of 1 that are determined curve-fitting the results of our simulation. 

We estimate the elastic force f s (t ) on the basis of hydrodynamic sound-wave theory [25]. In a normal fluid, sound 
waves propagate according to a relationship between the pressure and the velocity of the fluid. Let us denote a small 
change in the pressure from its equilibrium value by p', a typical velocity of the fluid particles in the wave by v, and the 
velocity of sound by c s . If the condition v <C c s is satisfied, we have a relationship p' = pc s v for a traveling plane wave, 
where p is the constant equilibrium density of the fluid. We assume that this relationship is also satisfied in fluidized 
granular media under the same condition v <C c s . The velocity of sound c s is on the order of the thermal velocity c. The 
density p is on the order of M/(Sc 2 /2g), where c 2 /2g is the first order estimation of the surface height of the granular 
gas as described above, and S represents the area of the base of the system. In the vicinity of the bottom plate, v may 
be approximated by the velocity of the bottom plate vo(f) = AoOcos(of). Since f s (t) corresponds to p'S at the bottom 
plate, we have 

f s (t) = GM-A (ocos((ot) 7 (6) 
c 

where a is a numerical factor on the order of 1 that is used as a curve-fit parameter when we compare our theoretical 
predictions with the results of simulations. The condition v <C c s is also written as vo <C c in the vicinity of the bottom 
plate. Hence, the maximum value of vo(f ), AqO), must be much small compared to c: AqCO <C c. Similar estimation of 
the pressure of sound wave has already been discussed in Ref. [8]. 
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As property of the random force, we assume stationary Gaussian white noise: 

(R(t))=0, (R(t)R(t'))=l8(t-t'). (7) 

We also assume that the intensity of the random force / is determined by the fluctuation-dissipation relation, which is 
expected to be satisfied when the system is close to equilibrium state: 

/ = 2MpLk B T. (8) 

We will discuss later how the fluctuation-dissipation relation is violated in the stationary state which deviates far from 
equilibrium. 

Substituting Eq. (4) into Eq. (3), we have the following linear Langevin equation for the COM: 

f = -rf(Z-Z) + (9, 

Analytical solution 

It is straightforward to obtain the analytical solution of the Langevin equation (9) of the form 

Z(t)-Z = A o £sin(ftW + 0)+ [' G{t-t')^Q-dt' + F ini {t), (10) 

J — oo IVl 



where 



and 

CO 2 

tan 9 = - - 



11(0 

respectively. The function G(t ) is given by 



(-|<0<o), (12) 



G{t) = ^-sin(coof), (13) 

(Oq 

where coq = (Q 2 - (ju/2) 2 ) 1 / 2 . The last term F ini (t) consists of those that depend on the initial conditions and 
vanish after a sufficient amount of time. Thus, the term is negligible when calculating long-time averages of physical 
quantities in the stationary state. 

NUMERICAL TEST OF THE LANGEVIN APPROACH 

We have performed event-driven MD simulations of one- and two-dimensional systems. The one-dimensional system 
consists of inelastic hard rods on a vibrating bottom plate. The two-dimensional system consists of inelastic hard 
disks on a thermal bottom wall. These two systems have been investigated in many previous numerical and theoretical 
works. In our two-dimensional simulations, we have used an efficient event-driven MD algorism developed by one of 
the authors [26]. Here we will demonstrate how our theory based on the Langevin equation can describe macroscopic 
properties of the granular fluids. 



ID system on a vibrating bottom plate 

We first study a one-dimensional granular fluid on a vibrating bottom plate [2, 3, 12]. We performed simulations 
systematically by changing the number of particles N (N = 10, 20, 100, 1000), the restitution coefficient r (r = 0.80 ~ 
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0.9999), and the maximum acceleration of the bottom plate T = A Q co 2 /g (T = 10 ~ 2560). The physical quantities are 
averaged over a sufficiently long period of time. When we compare our theory with simulation, we evaluate the global 
temperature T by kgT /2 = Ek, where Ek is the stationary value of the kinetic energy per particle defined as 



E K 



: lim 



M J0 



* 1 £ m . , 2 , 



(14) 



where v, is velocity of the i-th particle. Hence, the thermal velocity c is evaluated as c = yJkgT Jm = \J 2Ek Jin. 

Using the solution of the linear Langevin equation, we can calculate some macroscopic quantities. Let us first 
consider the amplitude £ of the oscillation of the COM. Theoretical prediction is given in Eq. (1 1). It must be noted 
here that we consider a fluidized state with the timescale (t/ T re /) 2 <C 1. This can be rewritten by substituting T = 2n/(0 
and x re i = (c / g) / jl as d) 2 3> 1, where a> is defined as a> = (Oc/g. In this limit, we expanded £ in Eq. (11) in terms of 



a)- 2 : 



CO 



(i + o (or 2 )). 



(15) 



Thus, it is predicted that the amplitude is inversely proportional to the scaled angular frequency of the vibration. 
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FIGURE 1. The amplitude £ of the oscillation of the center of mass as a function of the scaled angular frequency a) = coc/g 
of the bottom plate vibration. The accelerations F used in the simulations are T = 10,20,40,80, 160,320, and 640, except for the 
simulations that experienced an "inelastic collapse" (an infinite number of collisions in a finite time) [27, 28]. The dashed line 
corresponds to the theoretical prediction tr/ft) given in Eq.(15) with a curve-fit parameter a = 1.5. 

Figure 1 gives the simulation results for £ as a function of 6). We obtained good agreement between the theoretical 
prediction (15) and the simulation result when we chose the curve-fit parameter a = 1 .5. 

Secondly, we consider the power injection by the bottom plate Ph, which has been the subject of recent studies [4, 
8, 13, 14, 16]. It is defined by using the force exerted by the bottom plate and its velocity vq: 



P h = lim 



1 



tM 



F b (t)v (t)dt. 



t M ^oo f M Jq 

Substituting Fh using Eq. (3) into Eq. (16), and integrating by parts twice, we have 

l r'M ( d 2 z \ , l r'M 

P h = lim — / M—7T +Mg v (t)dt = -Mat 1 lim — / Z(t)v (t)dt. 

t M ^°° t M J0 V * / tM JO 



Then, substituting Eq. (10) into Eq. (17), we obtain 

dA co C0 2 ((0 2 -Q. 2 ) 



P b /MgA Q co = 



GAqCD 

2 c (oo 2 - D 2 ) 2 + (p> a) 2 ~ 2 ~ 



(l + 0(a>- 2 )). 



(16) 



(17) 



(18) 
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The results obtained by neglecting terms on the order of a>~ 2 coincide with the scaling predicted by kinetic theories [4, 
13, 14]: Pb ~ Mg (Aq(q) 2 jc. Figure 2 shows that the scaling relationship (18) with the same curve-fit parameter as 
Fig.l, a = 1.5, agrees well with the simulations in the region where Aq(o/c < 1. This result is consistent with the 
condition AqCd/c <C 1 required for the formula given by Eq. (6) to be valid. 
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FIGURE 2. The power injection by the bottom plate scaled by MgAoCO, as a function of c/AqCO. Here c is the thermal velocity 
evaluated as c = \/2EgJm. The accelerations T used in the simulations were T = 10,20,40,80, 160,320, and 640, except the 
simulations that experienced an inelastic collapse. The dashed line corresponds to the theoretical prediction (rr/2) (AqCo/c) given 
by Eq. (18) with a curve-fit parameter rr = 1.5. 

Thirdly, let us consider the power spectrum Iqm for the height of the COM. According to the Wiener-Khinchin 
theorem, this can be calculated analytically from the Fourier transform of the two-time correlation function Vcm(^) 
defined by 



WcM(t)= lim — / (SZ(t')SZ(t' +t))dt' , 

<M^~ tM JO 



(19) 



where 8Z(t) = Z(t) — Z and the brackets (• • • ) indicate an average over the random force R(t). Substituting Eq. (10) 
into Eq. (19) and performing the Fourier transform, we obtain 



c 2 ) 



((S(aV-<») + <5(&? + (B)) 



2/t 



(6 2 -&V ) 2 + (£aV) 2 ' 



(20) 



where dY is the angular frequency ©' scaled by g/c: &>' — co'c/g. Hence, our theory predicts that the power spectrum 
consists of two terms: the first term gives the delta-functional peak at the frequency of the bottom plate oscillation 
d) = (Oc/g. The second term represents a continuous spectrum. Indeed it has already been shown in Refs. [2, 3] that 
the power spectrum for the motion of the COM in a fiuidized state consists of a continuous spectrum and a sharp peak 
at the frequency of the vibration. In Fig. 3, we present the scaled power spectrum for the case of = 100, r=0.99, and 
r=160 obtained from our simulation. The result of simulation indeed shows the sharp peak corresponding to the first 
term and the continuous part corresponding to the second term. The part of continuous spectrum agrees well with our 
theoretical prediction when we choose the curve-fit parameters jx = 2.0 and Q. = 1 .5. 

The continuous part of the scaled power spectrum given by the second term of Eq. (20) predicts a universal behavior 
independent of the system parameters. In order to confirm this property, we first divided the logarithm of the angular 
frequency log ©' into bins with a constant interval, and took an average of the power spectrum over the bins; then 
we scaled it as given by Eq. (20). The results of simulations with various system parameters are shown in Fig.4. The 
simulation data indeed collapse on a single master curve which agrees very well with the curve of the theoretical 
prediction. 
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FIGURE 3. The power spectrum Iqm for the height of the center of mass scaled by c 5 /Ng^ obtained from simulations for 
N = 100, r=0.99, and T=160. The gray solid line depicts the theoretical prediction given by the second term in Eq. (20) with 
curve-fit parameters (L = 2.0 and O. = 1.5. 
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FIGURE 4. The power spectrum Icm f° r the height of the center of mass scaled by c 5 /Ng^ for various system parameters. The 
gray solid line is the same as in Fig. 3. 



Furthermore, we study how the law of equipartition fails when the inelasticity is increased. Let us define the 
stationary value of the kinetic energy of the COM for the motion along the z-axis as 



1 f'M \ 

T \ o 1 
tM^°° tM Jo \2 



K CM = lim — / ( -MV z (t) 2 )dt. 



Using Eq. (10), K C m can t> e calculated and yields 



M 
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(21) 



(22) 
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Figure 5 shows comparison of simulation results with the theoretical prediction (22). Here, the vertical axis shows 
the ratio of the kinetic energy of the COM, Kcm, to the kinetic energy per particle Ek- When the law of equipartition is 
satisfied, this ratio is equal to 1. The simulation data show deviation from the equipartition at Ek <C 10. This deviation 
can be well explained by our theory as far as it is not too large, using the same curve-fit parameter a — 1 .5 as used in 
Fig.l and Fig. 2. 
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FIGURE 5. The ratio of kinetic energy of the center of mass Kqm and kinetic energy per particle Equ for A' = 100 and different 
r values (from 10 to 2560), except the simulations that experienced an inelastic collapse. The gray solid line gives the theoretical 
prediction Eq. (22) with a curve-fit parameter a = 1.5. 



2D system on a thermal wall 

We then study a two-dimensional granular fluid on a thermal bottom wall [10, 11,21]. The system considered here 
consists of N inelastic hard disks with mass m and diameter d, in a two-dimensional space of a width L. We set periodic 
boundary condition in the horizontal direction and infinite height in the vertical direction. The stochastic thermal wall 
with the temperature 7b is set at the bottom. When a disk collides with the wall, it comes off with the value of the 
vertical velocity component, v„, sampled from the probability density 

KVn) = ^L exp f_^LV (23 ) 
1 ' k B T Q y \ 2k B T ) 

In order to prevent unphysical flows in the horizontal direction which may be caused by the periodic boundary 
conditions, we impose that the velocity component parallel to the wall remains unchanged by the collision with 
the wall. The dynamics is evolved by binary inelastic collisions between disks with a constant normal restitution 
coefficient r. The system can be completely characterized by the number of disks N, the dimensionless width L = L/d, 
the dimensionless driving intensity A = ksTo/mgd, and the restitution coefficient r [10, 11]. The value of L is chosen 
sufficiently small to prevent the density instability, so that the system remains homogeneous in the horizontal direction. 

Our theory can be easily applied to the case of the thermal bottom wall by simply modifying f s (t) — in Eq. (9). 
The other difference from the one-dimensional system is that the thermal velocity c is defined by c = \j2ksT Jm. 
When we compare our theory with simulation, we evaluate T by kgT = Ek and c by c = \J 2Ek Jm. 

Then, the theoretical curve of the power spectrum 7cm for the height of the COM is given by the second term in 
Eq. (20) divided by the system dimensionality 2, which comes from the change in the relation between kgT and Ek ■ 

Icm{co')/4-, - -- „Z ,, — • ( 24 ) 



N 8 (£2 2 - a' ) 2 + (A<» 



^2 
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In the Figs. 6 and 7, the power spectra Iqm of the COM with fixed parameters (L, A) = (10, 10 3 ) for various (N,r) 
are shown. Figure 6 shows the results of the simulations in the case of nearly elastic collisions (1 — r <C 1). Here we 
find that the scaled power spectra concentrate into a master curve. The theoretical curve with the curve-fit numerical 
factors £i = 2.0 and Q = 1.5 shows good agreement with numerical simulations. On the other hand, Fig. 7 shows the 
result of the simulations in the large inelasticity case; it demonstrates a systematic deviation from the curve predicted 
by our theory. We speculate that this deviation is caused by violation of the fluctuation-dissipation relation, which will 
be discussed in the following sections. 
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FIGURE 6. The power spectrum Iqm for the height of the center of mass scaled by c 5 /Afg 3 obtained from simulations for various 
(N,r). The driving intensity is A = 1000. Here the restitution coefficients r used in the simulations are close to 1 (nearly elastic 
case). The thick solid line depicts the theoretical prediction given by Eq. (24) with curve-fit parameters fl = 2.0 and £2 = 1.5. 




FIGURE 7. The same result as in Fig. 6 but the restitution coefficients r are smaller than those used in Fig. 6 (large inelasticity 
case). The thick solid line depicts the theoretical prediction given by Eq. (24) with curve-fit parameters fL = 2.0 and £1 = 1.5. 
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Failure of the law of equipartition in 2D system 



We first investigate how the law of equipartition fails when inelasticity is increased. Figure 8 shows the ratio of the 
kinetic energy of the COM K C m to ksT '/2(= E K /2). Since K CM is defined by Eq. (21) using only the z-component 
of the COM velocity V z , this ratio is to be equal to 1 when the law of equipartition is satisfied. We found a systematic 
deviation from the law of equipartition at low temperature; the deviation seems to be inversely proportional to T and 
independent of N. Our theory can not give any explanation on this behavior of the ratio of the energies. 



100 e 




0.01 0.1 



FIGURE 8. The ratio of kinetic energy of the center of mass Kqm and the half of kinetic energy per particle k B T /2 plotted against 
k B T/2. Here k B T is evaluated by k B T = E K . The restitution coefficients r are from 0.964 to 0.9995 for N = 1000, from 0.991 to 
0.9999 for N = 500, and from 0.90 to 0.999 for N = 200. The dashed line gives the law of equipartition of energy. The solid line 
gives a slope of — 1 .07. 

It is known from previous studies that the system shows density inversion in a certain condition [10, 11, 20, 21]. 
Here we discuss on a connection between the density inversion phenomena and the failure of the energy equipartition. 
In the previous study based on hydrodynamic equations, Bromberg et al. [21] have shown that if a parameter X defined 
by 

X^^N h (l-r 2 f\ (25) 

where Nf, is the number of layers of disks at rest, becomes larger than the critical value X c ~ 1.06569..., then density 
inversion develops in the system. Here we show in Fig. 9 the ratio of the energies as a function of X. We observe 
indeed that the deviation from the equipartition starts at a certain value of X near X c . We also confirmed that at the data 
points where the law of equipartition fails, the density inversion indeed develops as shown in Fig. 10. 



Failure of fluctuation-dissipation relation in 2D system 

The fluctuation-dissipation relation (of the second kind) (8) with respect to our Langevin equation (9) with f s = 
can be derived using the same procedure in the study of Brownian motion based on the Langevin equation [24]. To 
begin with, let us define the intensity of the random force / by Eq. (7). Then, we calculate K CM using the solution (10) 
of the Langevin equation (9) with f s (t) = 0, and represent it in terms of /: 

* CM= 4M^ (26) 

Finally, if the equipartition of energy, that is , Kqm = k B T /2, is satisfied, then we obtain the fluctuation-dissipation 
relation as given in Eq. (8). Therefore, in the case of large inelasticity where the law of equipartition of energy fails, 
the fluctuation-dissipation relation also fails to be satisfied. 
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FIGURE 9. The ratio of kinetic energy of the center of mass Kqm and the half of kinetic energy per particle kg T /2 plotted against 
X defined by Eq. (25). The dashed line gives the law of equipartition of energy. The dotted line indicates X c , the critical value of X 
for the density inversion. The solid line gives a slope of +2. 
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FIGURE 10. Snapshots and packing fraction profiles of simulations for A = 1000, N = 500, and different X values. 



Having observed that the fluctuation-dissipation relation is violated in the large inelasticity case, one can then go on 
to consider whether the linear Langevin equation itself is valid in the case of large inelasticity. In order to clarify this 
point, we use Eq. (26) and define / in terms of Kcm without use of the fluctuation-dissipation relation: / = 4M/j.Kcm- 
Then, the theoretical expression of I CM becomes 



W©')/ 



Kcm 
M 



(Q. 2 - co' 2 ) 2 + (A© 



(27) 
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When we compare Eq. (27) with the results of simulation, we determine Kcm from the simulation data. Figure 1 1 
shows the rescaling of the power spectra successfully collapses the simulation data onto a single curve near the peak. 
The curve of the theoretical expression (27) with fi = 0.6 and Q, = 1.5 seems to show better fitting than the case 
of fi — 2.0 and Cl = 1.5. This change in the coefficient ft in the case of large inelasticity may be attributed to the 
internal structure caused by the density inversion. We have not yet understand how the values of the phenomenological 
constants in our theory are related to the density profile of the system. 



100 1 — — — 




of 

FIGURE 11. The power spectrum 1cm f° r the height of the center of mass rescaled by 4(c / g) 3 Kcm /M for various (N,r). The 
parameters and lines of the simulations are the same as in Fig.7. The thick solid line gives the theoretical prediction given by 
Eq. (24) with curve-fit parameters fi = 2.0 and Cl = 1.5. The thick dashed line gives the theoretical prediction with fi = 0.6 and 
6= 1.5. 



CONCLUDING REMARKS 

A Langevin equation for the motion of the center of mass was formulated to describe macroscopic properties of 
the fluidized state of granular matter under gravity. The analytical expressions of some macroscopic quantities were 
derived and compared with the results of extensive simulations. There are three numerical factors, a, fx and Ci, that can 
not be specified by our phenomenological consideration; they were used as curve-fit parameters when we compared 
our theoretical predictions with the results of simulations. In the present study, we performed event-driven molecular 
dynamics simulations for the following two systems: a one-dimensional system of inelastic rods on a vibrating bottom 
plate and a two-dimensional system of inelastic disks on a thermal bottom wall. 

In the one-dimensional system, we found good agreement between our theory and simulation. The results of 
simulations suggest the following values of the numerical factors: a — 1.5, jX = 2.0 and Q. = 1.5. Deviation from the 
law of equipartition is observed when the kinetic energy per particle is sufficiently small. This deviation is explained 
well by our theory as far as it is not too large to violate the fluctuation-dissipation relation. 

In the two-dimensional system, the results of simulations are in good agreement with the theoretical predictions if 
the collisions between disks are nearly elastic. The results of simulations suggest the following values of the numerical 
factors: fi = 2.0 and Q. = 1.5. As the inelasticity of the interparticle collisions increases, however, the power spectrum 
for the center of mass obtained by the simulations gradually deviates from the prediction of theoretical curve. We 
speculated that this deviation is caused by violation of the fluctuation-dissipation relation. From the systematic event- 
driven simulation for the wide range of parameters, we found failure of the law of equipartition of energy with respect 
to the kinetic energy of the center of mass, which necessarily causes violation of the fluctuation-dissipation relation. 
Contrary to the one-dimensional system, the deviation from the law of equipartition can not be explained within the 
framework of our theory. Furthermore, we obtained evidence that this failure of the law of equipartition of energy 
is caused by the inversion of the density profile [10, 11, 20, 21]. Finally, we showed that the theoretical predictions 
obtained without assuming the fluctuation-dissipation relation agree with the results of simulations by choosing the 
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curve-fit parameters jl = 0.6, Q. = 1.5. This result suggests that the Langevin equation itself might remain valid even 
when the fluctuation-dissipation relation is violated, and it might be still useful to explain a universal behavior of the 
power spectrum in the large inelasticity case. 
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